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Abstract 

Factorization  of  linear  programming  (LP)  models  enables  a  large  portion  of  the  LP  tableau  to  be 
represented  implicitly  and  generated  from  the  remaining  explicit  part.  Dynamic  factorization  admits 
algebraic  elements  which  change  in  dimension  during  the  course  of  solution.  A  unifying  mathematical 
framework  for  dynamic  row  factorization  is  presented  with  three  algorithms  which  derive  from 
different  LP  model  row  structures:  generalized  upper  bound  rows,  pure  network  rows,  and  generalized 
network  rows.  Each  of  these  structures  is  a  generalization  of  its  predecessors,  and  each  corresponding 
algorithm  exhibits  just  enough  additional  richness  to  accommodate  the  structure  at  hand  within  the 
unified  framework.  Implementation  and  computational  results  are  presented  for  a  variety  of  real-world 
models.  These  results  suggest  that  each  of  these  algorithms  is  superior  to  the  traditional,  non-factorized 
approach,  with  the  degree  of  improvement  depending  upon  the  size  and  quality  of  the  row  factorization 
identified. 
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1.  Introduction 

A  recurring  theme  in  the  development  of  algorithms  for  linear  programming  has  been 
the  identification  and  exploitation  of  special  problem  structure.  Ideas  as  apparently  disparate 
as  the  bounded- variable  simplex  method,  primal  and  dual  decomposition  methods,  pure  and 
generalized  network  primal  simplex  algorithms,  primal  partitioning  and  column  generation 
schemes  may  be  unified  to  a  degree  with  this  view. 

The  factorization  approach  introduced  by  Graves  and  McBride  (1976)  isolates  special 
structure  in  LP  tableaus.  We  are  interested  in  using  factorization  to  reinterpret  existing 
algorithms,  and  to  discover  common  principles  and  apply  them  to  develop  new  algorithms. 
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Although  all  algorithms  developed  this  way  will,  in  theory,  solve  any  LP,  the  efficiency  of 
any  particular  factorization  approach  will  be  influenced  by  the  relative  number  of  factored 
constraints  and  their  influence  on  the  algorithm:  the  size  and  quality  of  the  special  structure 
isolated  determines  the  influence  of  any  particular  factorization  applied  to  any  particular 
LP. 

Based  on  prior  work  by  Brown  and  Graves  ( 1975),  in  which  generalized  upper  bound 
rows  were  successfully  incorporated  in  a  large-scale  optimization  system,  we  are  interested 
in  pursuing  dynamic  row  factorization,  where  the  dimension  of  the  factored  structure  may 
vary  (or  even  fail  to  be  present)  as  the  solution  progresses.  In  our  setting,  we  require  the 
row  structure  of  the  model  instance  to  be  specified  prior  to  solution,  and  that  this  structure 
remain  fixed  during  solution.  An  extension  of  this  approach  is  to  allow  the  row  structure  to 
vary  as  the  model  is  solved:  this  is  a  conceptually  simple  extension  of  the  approach. 

Each  algorithm  is  developed  by  factoring  the  constraints  of  the  LP  model  into  two  classes: 
those  that  have  the  special  structure  (factored )  and  those  that  do  not  (explicit).  This 
constraint  factorization  induces  a  factored  structure  in  the  LP  tableaus  which  is  exploited 
computationally.  We  demonstrate  the  dynamic  factorization  approach  for  three  special 
structures: 

-  generalized  upper  bound  rows; 

-  pure  network  rows;  and 

-  generalized  network  rows. 

We  implement  each  of  the  factorization  algorithms  by  integrating  it  within  the  X-System 
( Brown  and  Graves,  1975). 

While  the  terms  “partitioning”  and  “factorization”  are  frequently  used  interchangeably 
in  the  literature,  we  observe  a  distinction  between  the  two  approaches.  We  consider  parti¬ 
tioning  methods  to  be  based  on  special  structure  in  the  original  problem  instance,  which 
need  not  induce  special  structure  in  the  LP  tableau  —  in  fact,  the  method  need  not  be 
tableau-based.  In  contrast,  factorization  methods  are  based  on  special  structure  which  occurs 
in  bases  and  thus  in  the  basic  tableau.  Thus,  we  classify  dual  decomposition  (Dantzig  and 
Wolfe,  1960),  primal  decomposition  (Benders,  1962),  and  primal  partitioning  (Rosen, 
1964)  as  examples  of  partitioning  methods. 

Perhaps  the  earliest  example  of  what  we  consider  factorization  is  the  treatment  of  simple 
upper  bounds  by  Dantzig  ( 1954)  and  ( 1963)  and,  independently,  by  Charnes  and  Lemke 
(1954).  They  observe  that  it  is  more  efficient  to  enforce  the  “logical”  upper  bound 
constraints  with  logical  tests  within  the  algorithm  rather  than  treat  them  explicitly  along 
with  other  ‘  ‘structural’  ’  constraints.  While  not  originally  presented  in  the  context  of  a  formal 
tableau  factorization,  the  approach  is  easily  viewed  as  such. 

The  mutual  primal-dual  method  of  Graves  ( 1965)  focuses  attention  on  the  special  role 
of  nonnegativity  constraints  in  linear  programming.  A  clear  distinction  is  drawn  between 
the  computational  convenience  of  treating  nonnegativity  constraints  implicitly  rather  than 


G.G.  Brown,  M.P.  Olson  /  Mathematical  Programming  64  { 1994)  17-51 


19 


explicitly  and  the  unambiguous  mathematical  equivalence  of  all  problem  constraints,  struc¬ 
tural  or  nonnegativity.  Emphasizing  the  special  importance  of  inequality  constraints,  the 
approach  yields  an  elegant  theory  and,  as  we  will  see.  efficient  implementations.  We  view 
this  algorithm  as  the  first  formal  example  of  factorization. 

A  similar  primal-dual  algorithm  is  presented  by  Balinski  and  Gomory  ( 1965).  Related 
work,  in  which  efforts  are  made  to  exclude  slacks  from  the  product-form  representation  of 
the  primal  basis,  includes  that  of  Zoutendijk  ( 1970)  and  Powell  ( 1975). 

Dantzig  and  Van  Slyke  ( 1 967 )  extend  the  earlier  work  for  simple  upper  bounds  and  lend 
a  more  structured  treatment  to  generalized  upper  bounds  ( GUB ) .  In  a  problem  with  p  GUB 
constraints  and  m  structural  constraints,  their  approach  requires  a  working  basis  of  dimension 
m  +  1 ,  a  considerable  savings  when  p  is  large. 

Hartman  and  Lasdon  ( 1 972)  specialize  this  approach  to  the  multicommodity  capacitated 
transshipment  problem.  In  this  case,  the  structure  of  the  basic  pure  network  columns  intro¬ 
duces  additional  structure  into  the  working  basis,  allowing  further  simplifications  in  basis 
representation  and  update  techniques.  Helgason  and  Kennington  ( 1977)  develop  techniques 
for  representing  the  working  basis  in  product  form  and  provide  graphic  interpretation  of  the 
basis  updates.  Kennington  ( 1977)  reports  an  implementation  of  the  algorithm. 

McBride  ( 1972)  and  Graves  and  McBride  ( 1976)  formalize  and  generalize  the  factori¬ 
zation  approach.  They  view  factorization  as  a  unifying  framework  for  tableau-based  simplex 
specializations  and  illustrate  this  by  developing  a  variation  of  the  GUB  algorithm  of  Dantzig 
and  Van  Slyke  and  a  GUB  algorithm  for  the  doubly-coupled  linear  programs  of  Hartman 
and  Lasdon  (1970).  They  present  a  new  algorithm  for  the  set  partitioning  LP  and  an  equality- 
constrained  form  of  the  pure  network  with  side  constraints  model.  Brown  and  Graves  (1975) 
report  an  implementation  of  inequality-form,  dynamic  GUB  row  factorization  for  large- 
scale  problems. 

Schrage  ( 1975)  extends  the  succession  of  simple  and  generalized  upper  bounds  by 
introducing  variable  upper  bounds  ( VUB ) ,  which  are  constraints  of  the  form  Xj  <  xk,  where 
jc*  is  said  to  be  the  variable  upper  bound  of  Xj.  Schrage  implicitly  represents  the  VUB 
constraints  by  expressing  VUB  variables  in  terms  of  other  variables.  This  permits  the  basis 
representation  to  be  treated  in  two  parts,  one  a  large  matrix  which  changes  infrequently  and 
thus  needs  only  occasional  update,  and  the  other  a  small  working  basis  which  requires 
regular  attention.  Thus,  computation  and  storage  savings  may  be  realized.  Schrage  ( 1978) 
extends  these  ideas  to  what  he  calls  generalized  VUB  (GVUB)  constraints,  which  arise 
frequently  in  models  with  fixed  charges. 

Klingman  and  Russell  (  1975)  sketch  a  factorization  method  for  solving  transportation 
problems  with  side  constraints.  They  suggest  techniques  for  performing  simplex  iterations 
and  updating  the  problem  representation.  Chen  and  Saigal  ( 1977)  present  a  similar  approach 
for  solving  capacitated  network  flow  problems  with  additional  linear  constraints.  Both  of 
these  presentations  employ  a  graphical  description  of  the  basis  update  and  treat  the  basis  in 
two  parts:  one  corresponding  to  a  rooted  spanning  tree  defined  on  the  underlying  graph,  and 
the  other  a  general  working  basis.  Glover  et  al.  ( 1978)  report  an  implementation  of  the 
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Klingman  and  Russell  design,  but  one  which  (curiously)  only  accommodates  a  single  side 
constraint.  McBride  ( 1989)  reports  an  implementation  which  requires  the  pure  network 
rows  to  be  equalities  and  allows  more  than  one  side  constraint. 

Generalized  networks  with  side  constraints  are  addressed  by  Hultz  and  Klingman  (1976), 
who  present  details  for  the  simplex  priceout,  column  generation,  and  basis  update.  Hultz 
and  Klingman  (1978)  report  an  implementation  that  (curiously)  solves  the  “singularly 
constrained”  generalized  network  problem.  McBride  (1989)  reports  an  implementation 
that  is  not  restricted  to  a  single  side  constraint. 

The  factorization  approach  has  been  extended  by  consideration  of  embedded  structures. 
Glover  and  Klingman  ( 1981 )  consider  an  LP  with  embedded  pure  network  structure,  i.e., 
the  pure  network  structure  appears  in  only  a  subset  of  the  rows  and  columns  of  the  tech¬ 
nological  coefficient  matrix.  They  give  an  algorithm  similar  in  spirit  to  their  pure  network 
with  side  constraint  model,  but  the  presence  of  the  “side  variables”  significantly  compli¬ 
cates  the  basis  representation  and  update.  They  report  an  implementation  of  the  algorithm 
but  (curiously)  restrict  test  problems  to  have  no  complicating  variables. 

McBride  ( 1985 )  solves  an  LP  with  embedded  generalized  network  structure,  presenting 
methods  for  pricing,  column  generation,  basis  representation  update  and  data  structures.  A 
successful  implementation  is  reported  to  be  about  five  times  faster  than  MINOS  (ca.  1977: 
Murtagh  and  Saunders,  1977)  for  the  models  tested. 

Algorithms  to  solve  problems  with  special  substructures  have  motivated  research  to 
efficiently  identify  such  substructures.  Brearley,  Mitra  and  Williams  ( 1975)  describe  algo¬ 
rithms  for  detecting  GUB  row  sets  and  exclusive-row  structure  sets  (a  set  of  rows  whose 
structure  may  be  transformed  to  GUB  by  column  scaling).  Greenberg  and  Rarick  ( 1974) 
and  Brown  and  Thomen  ( 1980)  develop  algorithms  to  identify  GUB  sets.  Brown  and 
Wright  ( 1984)  identify  pure  network  constraint  substructures.  Brown,  McBride  and  Wood 
( 1985)  present  a  method  for  locating  embedded  and  row-only  generalized  network  struc¬ 
tures. 

Todd  ( 1983)  develops  a  geometric  interpretation  of  factorization  which  is  for  our  pur¬ 
poses  equivalent  to  the  algebraic  development  of  Graves  and  McBride  ( 1976). 

In  the  following  sections  we  establish  notational  conventions,  develop  the  mathematical 
foundation  of  a  primal-dual  simplex  method  and  show  the  effects  of  row-factorization. 
Next,  a  general  row-factorization  algorithm  is  developed,  specialized  to  GUB,  pure  network, 
and  generalized  network  rows,  and  tested  on  a  suite  of  real-world  problems.  Finally,  we 
discuss  how  the  methods  can  be  generalized  further  and  how  they  can  be  applied  more 
effectively. 


2.  Mathematical  preliminaries 


The  traditional  statement  of  the  linear  programming  (LP)  problem  is 
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(LP)  min  wy 


s.t.  i=l,...,  m, 

e}y> 0,  j=l,...,n, 

where  v  is  an  n-vector  of  decision  variables,  w  a  vector  of  cost  coefficients,  each  a,  an  rt- 
vector  of  technological  transformation  coefficients,  each  r,  a  scalar  right-hand  side  coeffi¬ 
cient,  and  ej  the  y'th  unit  vector.  While  this  statement  of  the  problem  is  clear  and  unambiguous, 
there  are  reasons  for  preferring  an  alternative.  The  insistence  upon  drawing  a  formal  dis¬ 
tinction  between  the  “structural”  constraints  a,v  <  r,  and  the  “nonnegativity”  constraints 
e}y  >  0  obscures  the  mathematical  structure  of  the  problem  by  suggesting  that  the  two  types 
of  constraints  are  inherently  different.  Certainly  the  exploitation  of  the  special  structure  of 
the  ejy  >  0  constraints  leads  to  computational  efficiencies;  however,  in  our  theoretical  devel¬ 
opment  of  the  algorithm,  we  prefer  to  treat  them  simply  as  general  inequality  constraints. 

In  order  to  achieve  a  consistent  form,  we  rewrite  the  nonnegativity  constraints  as  —  eyy  <  0 
and  group  them  with  the  structural  constraints.  The  problem  statement  then  becomes 

(PLP)  min  wy 

V 

s.t.  a,y<r;,  i~  1 . m  +  n , 


where  wy  is  called  the  extremal  function. 

From  the  standpoint  of  a  primal  algorithm,  a  matrix  partitioned  form  of  the  primal  tableau 
is  derived.  Let  {a„ ,  ai2, .. .,  ain }  be  a  basis  for  IR”  at  y°.  For  notational  convenience  we  will 
partition  the  constraints  into  two  sets,  those  that  are  basic  (binding)  aty°  and  those  that  are 
nonbasic  (not  necessarily  binding)  aty°. 


b\  " 

"  ah  " 

7." 

~rh  " 

B  = 

b2 

«/: 

.  /= 

fz 

ri2 

_ b 

Jn_ 

_n„  _ 

~d,  " 

Q in  + 1 

8 1  " 

r. 

1  hi  +  1 

D  = 

d2 

- 

^ in  +  2 

II 

&0 

8z 

- 

Y 

r  hi  +  2 

_dm  _ 

_“l’n+in  _ 

_§i)1  _ 

Y 

_  f«  +  m  _ 

Using  this  notation,  the  current  basic  solution  y°  may  be  expressed  as  By0  =f  and  since 
the  rows  of  B  are  by  definition  linearly  independent,  y°  exists. 

To  isolate  the  important  algebraic  components  let  us  assume  that  at  the  current  basic 
solution  the  basis  consists  of  h  structural  constraints  and  n  —  h  nonnegativity  constraints. 
Then 
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In  partitioned  matrix  form, 

h  n~h 

B  =  (au  Ai2  \h  , 

\  0  -I'n-h 

and  thus 

h  n  —  h 

p=-b~'=Ia;, -Ar,'A12^A 
l  0  /  Jn-h 


Similarly,  D  can  be  written  as 


and  in  partitioned  matrix  form 


h  n-h 


and  we  shall  call  DP  the  principal  part  of  the  tableau.  By  partitioning  w-  ( w,,  w2),gr  =  (#i, 
g:)T  and  y°  =  (y?,  y? ) ,  the  complete  tableau  may  be  written  in  partitioned  matrix  form  as 
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h 

n-h 

1 

A-' 

Au'Al2 

v° 

>  1 

h 

“A21A  n1 

~  A2)A  1 1 1  A|2 

S2  A  2  i_v  1  A22V2 

m  —  h 

— vv ,  a  n 1 

vv2 

—  vv,A  i~i 1 A ,  2 

O 

* 

1 

] 

Note  that  y?  is  displayed  explicitly  in  this  tableau.  Also,  y?  =0  since  the  corresponding 
nonnegativity  constraints  are  basic  and  thus  binding. 

The  corresponding  dual  problem  is 

(DLP)  max  .xr 


s.t.  x aj^wj,  j=\,...,n, 

xe,<0,  i=  1 . m. 

To  develop  a  matrix  partitioned  form  of  the  dual  tableau,  we  proceed  as  before.  Assuming 
the  dual  basis  consists  of  h  structural  constraints  and  m-h  nonnegativity  constraints,  we 
have 

T=  (/',  t2 . tm)  =  (aJ\  aj- . aj",  eiH*\  <?""), 

«  =  («',  M“,  «'")  =  O'1,  wh,  wjh,  0,  0), 


so 

H=(«‘,  u2)  =  (w',  0). 

The  nonbasic  constraints  are  then 


K=(k\ k2,  k")  =  (e‘\  eh- . e"\  aih+' . aJ"), 

r=(v\ t'2,  ....  f")  =  (0,  0,  ...,  0,  wJh",  ...,  h>), 


so  r=  (t’1,  v2)  =  (0,  w2). 

The  matrix  partitioned  form  of  the  basis  is  then 


h  m  —  h 

t=u ,,  o  y  , 

^  A21  /  7m— h 

and  with  the  choice  of  Q  -  T  ~ 1 , 

h  in  —  li 

a,v  0  y  , 

A2iAj"jl  /  / m-h 
with  the  remaining  constraints  forming 


t 
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The  principal  part  of  the  dual  tableau  is 


A  ~ 1 

A  U 

—  4;i  A  1 1 1 


0  /  Al2  _  A,,1  An'A^ 

/  0  A  22  _  A  2 1  A  j  j  ^  A  2  2  “  A  2 1  A  ||  A]2 


which  we  find  to  be  exactly  the  principal  part  of  the  primal  tableau,  so 


DP  =  QK. 


Thus  a  single  tableau  representation  supports  both  primal  and  dual  algorithms. 


3.  Interpretation  of  primal  and  dual  forms 

We  may  interpret  a  primal  or  dual  algorithm  as  simply  different  perspectives  of  this  same 
tableau,  wherein  a  primal  algorithm  basis  change  is  viewed  as  exchanging  primal  constraints 
and  a  dual  basis  change  exchanges  dual  constraints.  The  classical  Simplex  Method  may  then 
be  interpreted  as  solving  (PLP)  using  the  dual  perspective.  That  the  classical  Simplex 
Method  is  naturally  interpreted  as  a  dual  algorithm  comes  as  a  surprise  to  the  conventionally 
trained.  However,  the  consequent  mathematical  insight  is  compelling,  especially  in  light  of 
the  notational  simplification  and  apparent  underlying  role  of  A  I", 1 ,  which  we  refer  to  as  the 
transformation  kernel. 

There  are  several  reasons  for  preferring  a  primal-dual  algorithm  to  the  Simplex  Method. 
From  a  computational  standpoint,  because  slack  variables  are  carried  logically  rather  than 
introduced  explicitly,  we  are  able  to  clearly  identify  the  essential  information  needed  to 
execute  the  algorithm.  The  matrix  A  p, 1  plays  a  key  role  in  the  calculation  of  the  tableau, 
and  the  entire  tableau  can  be  constructed  from  A  1  and  original  problem  data.  Since  A  fj 1 
is  a  submatrix  of  the  inverse,  T  ~ ',  used  by  the  Simplex  Method,  it  is  smaller  and  requires 
fewer  arithmetic  operations  to  update  than  does  T~l. 

A  second  advantage  of  a  primal-dual  algorithm  lies  in  the  flexibility  it  offers  for  spe¬ 
cialization  to  particular  problem  classes  or  structures.  Indeed,  it  is  the  special  structure  and 
simplicity  of  the  nonnegativity  constraints  that  motivate  the  development  of  the  algorithm 
in  the  first  place.  It  is  frequently  the  case  that  other  special  structures  can  be  identified  in 
classes  of  (PLP).  Examples  of  such  structures  include  simple  upper  bounds,  generalized 
upper  bounds,  variable  upper  bounds,  pure  and  generalized  network  substructures,  etc.  Such 
structure  may  be  “static”  in  that  its  nature  and  dimension  remains  fixed  throughout  the 
solution  process,  or  the  structure  may  be  “dynamic”  in  which  case  its  precise  nature  and/ 
or  dimension  may  vary  as  the  problem  is  solved.  Some  special  structures  may  be  more 
strongly  characterized  by  their  column  structure  and  others  by  their  row  structure.  The 
primal-dual  perspective  leads  naturally  to  explanations  of  the  implications  of  virtually  any 
such  problem  structure  and  greatly  simplifies  the  implementation  of  such  a  specialization. 

When  an  LP  appears  as  a  subproblem  in  a  more  sophisticated  solution  setting  (for 
example,  in  a  mixed  integer  programming  problem  or  a  nonlinear  programming  problem ), 
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the  row/column  symmetry  of  a  primal-dual  algorithm  is  of  critical  importance  in  special¬ 
izing  the  solution  approach.  The  inherent  symmetry  of  such  an  algorithm  permits  easy 
adaptation  to  branch-and-bound  and  cutting-plane  approaches  to  mixed  integer  program¬ 
ming,  to  column  generation  settings,  as  well  as  to  primal  and  dual  decomposition  techniques. 

We  believe  the  reason  for  this  flexibility  offered  by  the  algorithm  lies  in  its  more  complete 
mathematical  foundation.  There  is  a  natural  consistency  that  arises  from  the  choice  of  a 
vector  space  having  the  same  dimension  as  the  problem  variables  that  is  lacking  in  other 
approaches.  A  natural  geometric  interpretation  of  the  solution  trajectory  follows  directly 
from  this  development.  Incidental  issues  such  as  finding  an  initial  basic  feasible  solution 
and  dealing  with  degeneracy  are  resolved  constructively  in  this  mathematical  framework 
(Graves,  1965).  Other  approaches  resort  to  unnecessarily  complicated  tangential  efforts. 

All  the  research  results  reported  here  can  be  developed,  with  some  effort,  in  the  framework 
of  the  classical  Simplex  Method.  However,  we  choose  to  present  these  results  in  the  manner 
of  their  development  —  the  mutual  primal-dual  view  presented  by  Graves. 


4.  Column  and  row  generation 


Rather  than  maintain  a  complete  tableau  DP,  now  consider  the  generation  of  just  column 
c  of  this  tableau.  Rewriting  in  a  manner  that  highlights  our  intentions,  and  labelling  row 
and  column  partitions  for  identification 


DP=  (i)  ( 
(ii)  V 


(i) 

[An'] 

—  A2\ [A  i  i  1  ] 


(ii) 

[An'Al2] 

a22  — a2|[a  i  i  1  a  j23 


)■ 


By  properly  sequencing  our  computations  we  will  exploit  the  fact  that  region  (ii)  of  a 
given  column  is  simply  a  linear  combination  of  terms  in  region  ( i)  of  the  same  column. 

Assume  we  want  to  place  the  current  representation  of  column  c  into  a  work  array  z, 
which  we  partition  as  zT=  (<:},  cl)  to  correspond  to  regions  (i)  and  (ii).  Expressed  in 
terms  of  the  transformation  kernel  A  fj 1 ,  we  compute  column  c  as: 
if  c  is  in  (j), 

Zt=[An'V, 


and  then 


—  A2,U,]; 
or,  if  c  is  in  (jj), 

~i  =  [A  1 1 1  (A12)‘  ] , 
and  then 


~2(A22)‘  A2,[c,]. 

Then  the  current  representation  of  column  c  is  available  in  cT  =  (;|,  zj  )• 
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The  computation  of  row  r  of  the  tableau  proceeds  in  a  similar  manner.  We  now  view  the 
principal  part  of  the  tableau  as 

'  (j)  <jj) 

DP  —  (i)  [A  T.1]  [An']Al2 

(ii)  [  A21 A  j  j  ]  Aoo  ~f*[  Ao,An  ]Ajt 

v  ~  ■■  *  " 

If  we  want  to  place  the  current  representation  of  row  r  in  a  work  array  z  partitioned 
conformably  with  (j)  and  (jj)  as  z=  (£3,  £4),  we  compute: 
if  row  ris  in  (i), 

S3  =  [Au']r, 
and  then 

24  =  [ ’3]  ^  12  ’ 

or,  if  row  r  is  in  (ii), 

S3  =  [(  —  A2i)rA  u  1  ], 

and 


£4  —  (^23),-  +  [S3]A|2, 

and  the  current  representation  of  row  r  is  available  in  z  —  (£3,  z4). 

We  see  that  in  each  case  calculations  proceed  by  first  using  a  representation  of  A  jj 1  to 
compute  a  portion  of  the  row  or  column  and  then  using  this  initial  computation  and  original 
problem  data  to  compute  the  remaining  part.  We  will  discover  that  our  specializations 
extend  this  approach  by  introducing  additional  tableau  partitions  which  allow  this  compu¬ 
tational  strategy  to  be  applied  on  a  larger  scale. 


5.  Transformation  kernel  update 

The  dynamic  behavior  of  A  jj1  is  important.  We  see  from  the  primal  row  basis  B  and 
nonbasic  rows  D  that  the  dimension  of  A jj 1  corresponds  to  the  number  of  basic  structural 
constraints,  or,  equivalently,  to  the  number  of  nonbasic  nonnegativity  constraints  (recall 
that  if  a  nonnegativity  constraint  is  nonbasic  and  thus  nonbinding,  the  corresponding  variable 
may  possibly  be  nonzero).  Recalling  that  our  primal  view  of  a  basis  exchange  is  as  an 
exchange  of  constraints  between  B  and  D ,  we  see  that  one  of  four  cases  may  occur  during 
a  pivot: 

•  A  structural  constraint  enters  the  basis  B  and  a  structural  constraint  leaves  the  basis  and 
enters  D.  Since  the  number  of  basic  structural  constraints  (and  the  number  of  nonbasic 
nonnegativity  constraints )  remains  unchanged,  the  dimension  of  A  fj 1  is  unchanged.  A  pivot 
of  this  type  involves  a  row  in  region  (j)  of  B  and  a  row  in  region  (ii)  of  D,  and  thus  it 
corresponds  to  a  pivot  coordinate  in  the  location  ( (ii),  (j) )  of  the  tableau  DP. 
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•  A  nonnegativity  constraint  enters  the  basis  and  a  nonnegativity  constraint  leaves  the 
basis.  Again,  the  dimension  of  A  fj 1  remains  unchanged.  Since  this  pivot  involves  a  row  in 
region  ( jj )  of  B  and  a  row  in  region  ( i )  of  D,  the  corresponding  tableau  DP  pivot  coordinate 
lies  in  ((i),  (jj)). 

•  A  structural  constraint  enters  the  basis  and  a  nonnegativity  constraint  leaves  the  basis, 
and  thus  the  number  of  basic  structural  constraints  (equivalently,  the  number  of  nonbasic 
nonnegativity  constraints)  increases  by  one.  The  dimension  of  A,-,1  is  increased  by  one. 
This  corresponds  to  a  pivot  coordinate  in  region  ( (ii),  (jj) )  of  the  tableau  DP. 

•  A  nonnegativity  constraint  enters  the  basis  and  a  structural  constraint  leaves  the  basis, 
and  thus  the  dimension  of  A  fj 1  decreases  by  one.  The  corresponding  pivot  coordinate  in 
DP  is  ( (i),  (j) ). 

We  see  that  we  may  exert  some  influence  on  the  behavior  of  the  dimension  of  A  fj 1  by 
our  strategy  for  selecting  target  exchanges  for  primal  and  dual  constraints  (i.e.,  our  pricing 
strategy)  and  through  our  tie-breaking  rules  for  choosing  pivot  row/column,  and  that  this 
dynamism  is  an  inherent  feature  of  an  effective  algorithm.  We  have  already  seen  the 
fundamental  importance  of  the  kernel  (A, , )  in  our  computations.  Thus,  a  successful  imple¬ 
mentation  must  manage  this  dynamic  behavior  efficiently  and  reliably. 


6.  Factorization 

The  row-factorized  problem  to  be  considered  is 


min  wy 

y 

C/7 

$ 

/A 

(explicit  constraints), 

Fy^b 

(factored  constraints). 

—  7y«0 

(nonnegativity  constraints), 

where  y  is  an  n- vector  of  decision  variables,  w  a  vector  of  cost  coefficients,  E  a  matrix  of 
constraint  coefficients  for  ‘  ‘explicit’  ’  constraints  with  right-hand  side  m- vector,  r.  Fa  matrix 
of  constraint  coefficients  for  “factored”  constraints  with  right-hand  side  p- vector  b,  and 
-  /  the  negative  of  the  identity  matrix.  In  this  general  development,  we  refer  to  the  F-type 
constraints  as  “factored”  only  to  distinguish  them  from  the  “explicit”  £-type  constraints, 
and  assume  nothing  about  their  structure.  Not  until  our  specializations  later  will  we  impose 
special  structure  on  F,  and  the  structures  we  will  consider  may  permit  the  representation  of 
the  F-type  constraints  without  the  inversion  of  a  matrix.  We  will  see  that  this  approach  is 
centered  around  handling  the  part  of  the  basis  corresponding  to  the  F-type  constraints 
explicitly  while  factoring  the  portion  of  the  basis  corresponding  to  the  F-type  constraints. 
The  notation  is  chosen  to  suggest  this  idea. 

Recall  that  a  basis  for  the  primal  algorithm  consists  of  n  linearly  independent  rows  from 
the  constraint  matrix  when  it  is  assumed  to  include  both  structural  (explicit  and  factored) 


28 


G.G.  Brown,  M.P.  Olson  /  Mathematical  Programming  64  ( 19941  17-51 


and  nonnegativity  constraints.  Assume  that  the  current  row  basis  consists  of  k  rows  from  E, 
I  rows  from  F  and  n  —  (k  +  l)  rows  from  - /.  Repeating  our  notation 

k  +  l  n-[k  +  l ) 

B  =  (A\\  A 12  'j  k  +  l  , 

\  0  —  /  /  n-{k  +  l) 

where  [A! ,  A12]  includes  all  basic  structural  rows,  from  both  E  and  F. 

We  will  ultimately  be  interested  in  isolating  the  effect  of  each  type  of  structural  constraint 
algebraically  in  the  factored  tableau,  and  thus  we  require  greater  resolution  in  our  factored 
basis.  Introducing  obvious  notation,  we  have 

k  l  n-(k  +  l ) 

[An  A12]  — / Eu  E\2  El3  \k, 

v  F2 1  F22  F 23  /  / 

where  the  kernel  of  dimension  k+ /  is  given  by 


Because  A , ,  is  a  basis  for  R*  + '  it  follows  that  it  is  always  possible  to  identify  among  the 
columns  of  [F2]  F22]  a  nonsingular  submatrix  F22  of  dimension  /,  since  otherwise  the  rank 
of  [F2l  F22]  is  at  most  /—  1  and  thus  the  rank  of  A, ,  is  at  most£  +  /—  1,  or  equivalently  the 
rank  of  B  is  at  most  n—  1.  We  will  later  see  that  one  of  the  important  implementation 
challenges  is  the  task  of  efficiently  managing  the  structure  and  nonsingularity  of  F22- 
The  full  factored  row  basis  is  then 

k  I  n  —  (k  +  l ) 

B=  ,j)[£l1  £|2  £|'  V 

( jj )  F2\  F 22  F23  I 

(iij)  0  0  —I  n-(k  +  l) 

Introducing  the  notation 

Au=£n  E\2F  it1  F 2 1 ,  A 1 3  =E\  3  Ef  2  F  22 1 F23 , 

where  A, ,  is  the  Schur  complement,  or  Gauss  transform,  of  F22  in  An  (e.g.,  Golub  and  Van 
Loan,  1983),  we  can  write  its  inverse  as 

Am'  ~A\\XEX2F22  An'A^ 

B~'=  -F^F^A;,'  U+Fx'F2tAu'El2)F,V  F22'{F23-F2]Au'A]3)  . 

0  0  -/ 

Grouping  the  coefficients  of  the  nonbasic  constraints  and  applying  the  same  column 
ordering  yields 
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(i) 

D=  (ii) 

(iii) 

(iv) 


k  I 

-I  0 
0  -/ 
E3\  E)  2 
Ea]  Fa2 


n~  (k  +  l) 
0 
0 
£33 
Fa, 


I 

m  -  k 
I  p-l 


The  principal  part  of  the  factored  tableau  is  DP,  where  P=—B  1  is  the  conjugate  row 
basis.  With  the  additional  notation 


A3 !  — F3  ] 

~E32F2: i 

'  F2  i  , 

A33  — £33 

—  E32F  22 

F2  3, 

> 

III 

~FaiF  22 1 

f21, 

F 43  —F 43 

Fa2  F  22 ' 

F23, 

the  principal  part  of  the  factored  tableau  is 


(j) 

(ii) 

(iii) 

(•) 

4-1 

—  A\\E\2F  22 

DP=  (ii) 

~F 22 1 F 2IA  i  |' 

U+Fz2xF2iA^EX2)Ff2' 

F  22X(F23— F2\A\f  AX3) 

(iii) 

—  A31 A  | , 1 

~E32)F 22 

A33  —  A31A  n  1  Al3 

(iv) 

\  ~  Fa\  A  j  1 1 

F43  ~  F 4|A||IA]3  / 

Partitioning  w=  (wl5  w2,  w3,  w>4),  rT=  (r|,  rT)  and  bT  =  (bj,  bj)  and  introducing  the 
notation 


W2=W2 ~WiF2VF2I,  VV-,  =w3  —  w2  F 12 1 F23 , 

b2  =b2  F ^  \ F 22  b | ,  r i  =r x  E\\F22b\,  r2  =r2  E2\ F 22  bx , 

the  complete  factored  tableau  is 


r  a-' 

n  1 1 

-Au'E^FCA 

A7,'A[3 

z  1 

:?! 

_ 1 

-FZ'FuAu' 

V+F^F2iAHXEn)FZ} 

F  22 1  ( F 2?  ""  ^2 1 A  j  1 1 A 1 3 ) 

F 22'  (h\  ~  F2\A\  i*  n ) 

-As,Ar,' 

( Aii A\ j 1  E\2  ~ E32)F 22 1 

^33  ”  ^ 3 1  ^  1 1*^13 

Fi  ~ A3] A n’ T| 1 

-FA'n' 

(  F41A  n'  £|1  —  F 42  )F  22 1 

F 43  ~  F 41 A ]  ] 1  Aj3 

b2-FMA„'rx 

—  VV’jA  1 , 1 

<VrsA,V£|2  H1 1  )  £  22 * 

m’j  wsA  1 1  ^  A  |  3 

H’,  FcAb,  +  w2A  j", 1  f| 

We  see  that  with  knowledge  of  the  current  factorization,  we  can  construct  the  entire 
tableau  from  F  j;1 ,  A 1  and  the  original  problem  data.  The  dimension  of  F  22  is  equal  to 
the  number  of  F-type  constraints  that  are  currently  basic,  and  thus  can  be  at  most  p.  The 
dimension  of  A  |j 1  is  equal  to  the  number  of  F-type  constraints  that  are  currently  basic,  and 
thus  cannot  exceed  m.  We  call  A  M 1  the  explicit  transformation  kernel  and  F  f,'  the  factored 
transformation  kernel. 
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7.  Factored  column  and  row  generation 

Consider  generation  of  column  c  from  the  principal  part  of  the  tableau  DP.  Rewriting  in 
a  manner  that  highlights  our  intentions 

(i) 

DP=  (ii) 

(iii) 

(iv> 

where  A , 3  is  defined  as  before,  and  the  brackets  “  [  ] "  and  “  {  } '  ’  contain  terms  common 
to  but  displayed  only  once  for  each  column. 

Assume  we  want  to  place  column  c  into  work  array  z.  We  partition  z  conformably  as 
zT=  (;},  zT>  zj,  zj ),  refer  similarly  to  components  of  unit  vectored  and  employ  a  (m+p)- 
vector  work  array  z.  The  notation  “  denotes  simple  assignment,  “  =  ”  indicates  that  a 
set  of  factored  equations  must  be  solved,  and  “=>”  explains  the  corresponding  result. 

If  column  c  is  in  ( j ) , 

Zx<-[Au'V, 

and  then 

Z* 

£ 22<-2  =  Z  =*  ^2  (  —  F  22*  £j|  [Cl  ]  }> 

then 

~3  *  £3  j  [ £|  ]  —  Eyi  { C2  }  ■> 

and  finally 

-4  *  F4l[c,J  —  F 42(^2 }» 

if  c  is  in  ( jj),  solve 

F22Z2  =  -(ezy  =>  ;>  * - (F 22  )', 

Z*~Ei2Z2  =*  Z*  (  El2F  22*  )*  1 

c  1  *  ^4 1 1 1  c  =*  tii  [  —  A 1 1 1  ElzF 22']% 

then 

z*~  — ^2i[^i]' 

F22z2=z  =>  z2^{(F22'Y-F;1'F21[z,]}' 
and  finally 


{~F^F2l[  ]} 
]~E,2{  } 
-£4,1  l-£42{  } 


<JJ)  (JJJ> 

[  ~A ,/ E^Fii1] 

[F22-F21F2A  ]}  { £  22 '  £23  F  22  F 2.[  ]} 
~E3l[  ]-£, 2{  }  £33  £3 1  [  ]-£»{  } 

-£4 ,[  3  £42 {  3  £43  - £41  [  3  £42 {  3 
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Cl  £3  ]  [  C  i  ]  £*32  {  ) » 

c4  4 — £4 1  [  c:  1  ]  —  £42  { ~2 )  i 

if  c  is  in  { jjj ) , 
c4-  (£3)*. 

Fn  2  C;  =  C  =>  Z2  F  22 '  (  F22  )  ‘ , 

“  ^  (£53)'  £  1 2  -^2  =>  C4- (A]3)% 

Ci  ■*— A , , 1  c  =*  “  1 4  [ A 1 1 1 A , 3  ]  ‘ , 

then 

z*~  ( £ 23 )  £ 2 1  c ( - 

F22z2=z  =»  c2 { £ ^ 1  ( £2.3 ) ‘  — •  £ 22 1  £2 1  [ c  1  ] } » 
and  finally 

C3  ( £33 )  *  £3 1  [  c  i  ]  £32  {  C^2  }  •> 
c4  (F42)1  ~  F41  [Ci  ]  —  E42  {c2  }• 

Similarly,  row  r  can  be  generated.  The  principal  part  of  the  tableau  is  now  viewed  as 


(j) 

(ii) 

(jjj) 

[  ]£l3  +  {  )  f22 

(i) 

[Au'1 

{-[  ]£,2£2V} 

DP=  (ii) 

[ 

~  F  22*  F21A  |~ | 1  ] 

1  [F^-l  ]El2F£) 

[  ]£,  3  +  {  }F2  3 

(iii) 

[  _  A,  |  A  u  1  ] 

{ (  - £32  -  [  ]£,2)£2aI 

}  £33 +  [  ]£,3  +  {  )F23 

(iv) 

k 

[-F^AT/] 

{ {  -  £42  -  [  \El2)F£ 

}  £43  +  [  ]£,.3  +  {  )Fj 

where  A3 ,  and  £4 ,  are  defined  as  before,  and  the  brackets  “  [  ] "  and  “  {  }  ”  contain  terms 
common  to  but  displayed  only  once  for  each  row. 

Assume  we  want  to  place  row  r  into  work  array  ~.  We  partition  z  conformably  as  z —  (c5, 
C7 ) ,  refer  similarly  to  components  of  unit  vector  er ,  and  employ  a  ( m  +  p)-ve ctor  work 
array 

If  row  r  is  in  (i), 

C5  [An1], 

then 


En 


*•5  ^  1 2 


Cfi  F  22  c 


[A,, '],.£, 2, 

4~  {  l~i 1  J r^l2 £ 22*  1* 


and  finally 
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Z*"~  [Oi]£l3  +  {Zfi  }^23' 

if  r  is  in  (ii),  solve 

^6^22  ~~(eb)r  **  "6  *  (^22 ')/-* 

<■*_  -6^21  Z«  (^*22*  )/-^21  * 

^  -^II'  =*  Z.<s  «“  [  —  F  22 '  ^*2 1  ^  11  *  ]  r » 

then 

z*-(eb)r-  [c5]£i2, 

<■6^22  ~Z  ^  Zf,  f  22*  )r  —  [ C5 ] £  1 2 ^ 22 '  }■ 

and  finally 

Z7  [Ss]£|3  "*■  {^6  23  * 

if  r  is  in  (iii), 

z*-~(E32)r, 

Z(,F 22  =  Z  ^6  *  (  ^^2)  22*  > 

Z«-(£3l)r+4^2l  =»  2<-(A3l)r. 

Z5  ^~  [ z  A !  1 1  ]  =>  Z5  «“  [  —  A31 A 1 1 1  ]r, 

then 

Z  4  ( Eyy ) r  ~  [Z.<;]£i2» 

Zfi  ^22  =  Z  =»  Zf,  *“  {  (  —  (£32)/-  —  C  ZZ5  ]  £" 1 2  )  /^  22 1  h 

and  finally 

z7  <-(£33),  +  [z's]£i3  +  {Zfi}^; 

if  r  is  in  (iv), 

z< —  (^42)^. 

C6F22=Z  =>  Z6  < - (  /^42  )  22 '  ' 

Z*“  (^41  )r  ~Zf,F2\  =»  Z4-  (^4l)r« 

Zj^tz'Ar,1]  -  ZsM-F-nA'r.'Jr, 

Z*  (  ^42  )  r  [  Z3  ]  £"  1 2  ’ 

Z6/r22=Z  ^  Zfi  *-  {  (  ~  (  F42  )  r  ~  [  Z.s  ]  E 1 2  )  f  22 '  } ' 


then 


G.G.  Brown.  M.P.  Olson/ Mathematical  Programming  64  ( 1994)  17-51 


33 


and  finally 


7-1  (F43),.  +  [Cs]£|3  +  {^6  } ^23 ■ 


8.  The  complete  algorithm 

The  complete  algorithm  is  described  in  terms  of  abstract  functions  which  operate  on 
fundamental  data  structures. 

Tableau  management  requires  two  index  maps:  one  yielding  the  intrinsic  coordinate  in 
the  principal  part  of  the  tableau  for  each  original,  extrinsic  problem  row  or  column,  and  the 
other  its  inverse  map.  Intrinsic  arguments  are  shown  in  lower  case,  and  extrinsic  in  upper 
case.  Index_Exchange(  index  l,index2)  updates  these  maps  for  the  exchange  of  a  pair  of 
tableau  coordinates  index  1  and  index2. 

The  tableau  regions  are  successive  partitions  of  indices: 


Tableau 

region 

ENDING 

INDEX 

Contents  of  region 

(i) 

MEC 

basic  Columns  solving  Explicit  rows 

(ii) 

MFC 

basic  Columns  solving  Factored  rows 

(iii) 

MER 

nonbasic  Explicit  Rows 

(iv) 

m 

nonbasic  Factored  Rows 

(j) 

NER 

basic  Explicit  Rows 

(jj) 

NFR 

basic  Factored  Rows 

(jjj) 

m  +  n 

nonbasic  Columns 

Increment!  endingindex)  and  Decrement ( ending  index)  are  functions  to  modify 
these  ending  indices. 

GenerateRow!  row )  and  Generate_Column(  column)  place  numeric  values  of  a  tab¬ 
leau  row  or  column  in  ROWCOL(  ) ,  which  is  commensurate  with  the  tableau  dimensions. 

Using  ROWCOL(  ),  UpdateJRim(row,coI)  maintains  current  numeric  values  of  the 
right-hand  side  and  bottom  row  of  the  tableau  in  RIM(  ). 

The  explicit  transformation  kernel,  A  ,"j '  is  operated  on  by  functions  using  ROWCOL(  ) : 

Add  E-Kernel  Row  ( ROW ) , 

Add  E-Kernel  Column  ( COLUMN ) , 

Delete_E-KerneI_Row(  ROW) , 

Delete  E-Kernel  Column!  COLUMN ) , 

Replace  E-Kernel  Row(  REPLACED  ROW, REPLACING  ROW ) ,  and 
UpdateExplicitTransformationKernel,  the  pivotal  update. 


Table  8. 1 

Secondary  and  tertiary  tableau  exchanges 

(j) 

Delete_E-Kernel_Column(  ROW ) 
IndexExchange!  col.NER) 

Index  Exchange!  NER.NFR ) 

Decrement!  NER ) 

(i)  Decrement! NFRl: 

Delete  E-Kernel  Row!  COL ) 
Index_Exchange(  row.MEC) 

Index  Exchange!  MEC.MFC) 

Decrement!  MEC) 

Decrement!  MFC ) 

DeletcE-KernelColumn!  ROW ) 

Index  Exchange!  col.NER  > 

IndexExchange!  NER.NFR ) 

Decrement!  NER) 

Decrement!  NFR); 

( ii )  Find_E-Kernel_Column_for_Kcy ( ROW, KIN ) 
Delete  E-Kernel_Row<  KIN ) 

Index  Exchange!  kin.MEC) 

IndexExchanget  row.MFC ) 

Decrement!  MEC ) 

Decrement!  MFC ) 


(jj) 

Index  Exchanget  col, NFR ) 

Decrement!  NFR ) ; 

Find_F-Kerncl_Column  to  Remove!  COL.KOU ) 
GenerateRow!  kou) 

Replace  E-Kemel_Row(  ROW, KOU ) 
IndexExchanget  row.kou ) 

Index  Exchanget  row.MFC ) 

IndexExchanget  MFC.MER ) 

Decrement!  MFC) 

Decrement  ( M  ER ) 

Index  Exchange!  col  .NFR ) 

Decrement!  NFR ) ; 

IF  (Factored_Kernel  Singular! COL, ROW) )  THEN 
Find_F-Kernel_Column_to_Remove(  COL.KOU ) 
Find_E-Kernel_Column_for_Key(  ROW, KIN ) 
Generate  Row(kou) 

ReplaceE-KernelRow!  KIN.KOU ) 
Index_Exchangc(  kin, kou ) 

ENDIF 

Index  Exchanget  row.MFC ) 

Index  Exchanget  MFC.MER ) 

Decrement!  MFC ) 

Decrement!  MER ) 


U) 


(jil  > 


( no  more  exchanges ) 


t  no  more  column  exchanges ) ; 

IF  tFactored  Kernel  Singular!  ROW, COL) )  THEN 
Find_E-Kernel_Column_for_Key(  ROW.KIN ) 
ReplaceE-KernelRow!  KIN.ROW ) 

Index  Exchanget  kin. row) 

ENDIF 
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(iii)  (  no  more  exchanges) 


Delete_E-Kernel_Column(  ROW) 
IndexExchange!  col.NER ) 

Decrement!  NER ) : 

Delete_  E-Kernel_Row  ( COL ) 

( iv )  Find  E-Kernel  Column  for  Key!  ROW, KIN ) 
IndexExchange!  kin.MEC ) 

Decrement!  MEC) 

Increment!  MER) 

Index  Exchange!  row.MER ) 


Add_E-Kernel_Column(  COL) 

Increment!  NER); 

Index_Exchange(  col.NER ) 

Find_F-Kernel_  Column_to_Remove(  COL.KOU ) 
Generate_Row(  kou ) 

Add  E-Kernel  Row(  KOU ) 

Index  Exchange!  kou.MEC ) 

Increment!  MEC ) 

Index  Exchange!  row.MER) 

Decrement!  MER ) 

( no  more  column  exchanges); 

IF  (Factored  Kernel  Singular!  COL.ROW) )  THEN 
Find  F-Kernel_Column  to_Remove<  COL.KOU ) 
Generate  Row(  kou ) 

Find_E-Kernel_Column_for_Key  ( ROW, KIN ) 
Replace  E- Kernel  Row(  KIN.KOU) 
IndexExchange!  kin,kou ) 

ENDIF 


Add  E-Kernel_Column(  COL ) 

Increment!  NER ) 

Increment!  NFR) 

IndexExchange!  col.NFR ) 

Index  Exchange!  NFR.NER ) ; 

Add  E-Kernel  Row!  ROW ) 

Increment!  MEC ) 

Increment!  MFC) 

Index_Exchange(  row. MEC ) 

IndexExchange!  MEC.MFC ) 

Increment!  NFR) 

Index_Exchange(  col.NFR ) ; 

IF  (Factored  Kernel  Singular! COL.ROW))  THEN 
Find_E-Kcrnel_Column_for_Key  ( ROW,  KIN ) 
Replace_E-Kernel_Row(  KIN.ROW ) 
Index_Exchange(  kin.row ) 

ENDIF 

Increment!  MFC) 

Increment!  MER ) 

IndexExchange!  row.MER ) 

Index  Exchange!  MER.MFC) 


A  pivot  in  tableau  row  “row”  and  column  “col”  is  associated  with  problem  row  or  column  indices  ROW  and  COL,  and  begins  with  Generate_Row(row), 
Generatc_Column(col),  Update_Rim(row,col),  and  primary  Index  Exchange!  row.col).  Next,  depending  upon  the  tableau  region,  secondary  and  even  tertiary 
exchanges  shown  above  may  be  required  to  preserve  the  factorization  and  the  tableau  structure  so  the  pivot  can  then  be  completed  with  Update- 
_Factored_Kernel(  ROW.COL)  and  Update_Explicit_Transformation_Kernel. 


OJ 
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A,', 1  can  be  represented  in  any  way  that  suits  the  implementer  and  efficiently  supports 
these  functions.  Generally,  we  find  A  |"j 1  to  be  relatively  sparse,  and  report  here  results  using 
an  array  with  an  entry-point  for  each  inverse  row  and  column  which  accesses  a  stack  of 
orthogonally-linked  nonzero  inverse  elements.  We  have  also  implemented  a  dense  vector- 
processor  version,  and  the  design  issues  for  LU-based  schemes  are  given  by  Olson  ( 1 989 ) . 

The  factored  kernel,  F22,  is  manipulated  with: 

Factored_Kemel_Singular(ROW,COLUMN),  a  Boolean  function, 

Find  E-Kernel  Column  for  Key(ROW.COLUMN),  a  column  from  region  (i), 

Find_F-Kernel_Column_toJRemove(  ROW  .COLUMN),  a  column  from  region  (ii), 
and 

Update_Factored_Kernel( ROW, COLUMN),  the  pivotal  update. 

The  complete  abstract  algorithm  is: 

Step  0.  Initialize. 

Step  1.  Select  primal  or  dual  algorithm. 

Step  2.  Select  a  primal  (col)  or  dual  (row)  violation; 

STOP  if  current  solution  is  terminal,  or 
Generate  Column( col)  or  Generate  Row( row). 

Step  3.  By  ratio  test  with  RIM(  )  and  ROWCOL(  ), 
select  a  (row)  or  (col)  pivot  coordinate; 

STOP  if  current  solution  is  terminal,  or 
GenerateRow(row)  or  Generate_CoIumn(col). 

Step  4.  Update  Rim(row,col), 

primary  Index_Exchange(row,col), 

perform  secondary  and  tertiary  exchanges  (Table  8.1 ), 

Update  Factored  Kernel  ( ROW, COL ) , 

Update  Explicit  Transformation  Kernel, 

Go  to  Step  1. 

Functions  for  maintaining  the  factored  kernel  vary  with  the  factorization,  and  a  good 
implementation  will  exploit  these  differences.  However,  a  general  specification  will  suffice 
for  all  factorizations  here. 

We  are  exclusively  interested  in  performing  two  fundamental  operations: 

1 .  Solving  factored  equations  of  the  form 

F 22^2  t>2 


and 
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where  ~2  and  z2  are  unknown  and  b2  and  b2  are  rational  (not  necessarily  integer). 

2.  Restoring  F22  to  the  desired  form  of  F22  which  makes  the  factored  equations  above 
easy  to  solve,  where  F22  results  from  inflicting  F12  with  a  column  exchange,  a  column  and 
row  deletion,  a  column  and  row  addition,  or  a  row  exchange. 

Factored  Kernel  Singular!  LABEL  1.LABEL2)  predicts  whether  F22  will  be  singular 
if: 

1.  LABEL  1  gives  a  basic  column  in  a  basic  factored  ROW  for  which  COL  =  LABEL2 
would  be  exchanged; 

2.  LABEL  1  is  a  basic  factored  ROW  solved  by  basic  column  COL  =  LABEL2,  both  of 
which  would  be  removed  from  the  factored  kernel; 

3.  LABEL  1  is  a  nonbasic  factored  ROW,  which  would  be  added  to  the  factored  kernel 
with  column  COL=LABEL2;  or 

4.  LABEL  1  is  a  basic  factored  ROW  which  would  be  replaced  by  some  other  row 
LABEL2. 

The  first  case,  a  column  exchange,  is  equivalent  to  asking  whether  solving  F22z2  =  b2  with 
b2  equal  to  region  (ii)  of  the  proposed  entering  column  COL,  yields  a  nonzero  term 
associated  with  the  basic  factored  row  in  ( row ) .  This  is  easy  to  answer  for  the  factoriza¬ 
tions  we  discuss.  For  instance,  Brown  and  McBride  ( 1984)  show  that  back-solving  a 
“nearly-triangulated  generalized  network’’  basis  F22  only  as  far  as  ROW  will  suffice,  and 
that  this  can  be  implemented  as  a  traversal  of  the  “backpaths”  of  the  (zero,  one,  or  two) 
coefficients  in  b2  for  the  column  now  basic  in  ROW.  If  numerical  precision  is  an  issue, 
’2  ( row )  can  be  computed  during  this  search  and  tested  for  significance. 

The  second  case,  a  row  and  column  deletion,  is  trivial  because  the  resulting  F22  will 
always  be  nonsingular. 

The  third  case,  a  row  and  column  addition,  can  be  answered  by  adding  ROW  and  COL 
to  P'22  without  disturbing  its  desirable  special  structure,  thus  creating  an  instance  of  case 
one.  For  instance,  placing  ROW  first,  and  COL  last  in  F22  suffices  for  all  factorizations 
discussed  here. 

The  fourth  case,  a  row  exchange,  can  be  answered  by  applying  the  second  case,  deleting 
ROW  and  its  basic  column,  and  then  (perhaps  repeatedly)  applying  the  third  case,  adding 
the  row  with  index  COL  and  any  column  which  would  yield  a  nonsingular  F22. 

Find  F-Kernel  Column  to  Remove(ROW,COL)  given  basic  factored  row  ROW, 
returns  its  basic,  or  “key’’  column  COL. 

Find_E-Kernel_Column_for_Key(LABEL,COL),  given  either  a  nonbasic  factored 
ROW  =  LABEL,  or  a  basic  column  LABEL  solving  a  basic  factored  ROW,  searches  for  an 
acceptable  basic  column  COL  in  Explicit  Rows  (region  (i))  using  Factored- 
Kernel_Singular(  ROW,COL) . 

Update  Factored  Kernel!  LABEL  LLABEL2)  restores  F22  to  the  desired  form  of  F22, 
with  a  possible  increase  or  decrease  in  dimension.  Following  the  case-by-case  scheme  of 
Factored_Kernel_Singular.  we  pre-  and  post-process  the  factored  basis  representation  to 
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permit  use  of  a  single,  static  factorization  update  of  conventional  design.  Olson  ( 1989) 
pursues  this  in  considerable  detail. 

Table  8.2  displays  supporting  data  structures  for  these  functions. 


9.  Computational  experience 

The  factorization  methods  introduced  here  have  been  implemented  and  used  to  solve  a 
variety  of  existing  models  provided  by  our  colleagues.  With  their  help,  we  have  extracted 
suggested  problem  instances  from  a  diversity  of  decision  support  systems  on  host  computers 
ranging  from  mainframes  to  micro-computers.  Because  our  goal  is  to  test  factorization 
technology  in  isolation,  the  results  reported  here  are  achieved  without  benefit  of  any  model- 
specific  knowledge  or  tuning. 

However,  we  also  seek  to  develop  effective  modeling  tools  for  customized  use  in  devel¬ 
oping  and  refining  new  models.  To  this  end,  we  have  greatly  benefited  from  the  experience 
and  advice  of  our  colleagues,  and  we  include  some  discussion  of  modeler  guidance  and 
insight  along  with  the  numerical  results. 

Each  model  is  introduced  below  by  a  short  synopsis.  Multiple  instances  of  some  models 


Table  8.2 

Factorization  algorithm  data  structures 


Factorizations 

Data  Structure 

Size 

Use 

GUB,PN,GN 

RIM(  ) 

m  +  n 

current  tableau  right-hand  side,  bottom  row 

ROWCOL(  ) 

m  +  n 

current  tableau  row  and  column 

MSKRC(  ) 

m  +  n 

logical  mask  true  for  corresponding  nonzero  in 

ROWCOL  (  ) ,  false  otherwise 

LQRC(  ) 

m  +  n 

LIFO  queues  of  nonzero  row  and  column  coordinates  in 
ROWCOL  (  ) 

KEY(  ) 

m  +  n 

basic  column  in  basic  factored  row,  and  vice  versa 

PN,GN 

WORK(  ) 

P 

values  for  bz  or  b2  in  basic  factored  equations 

MSKWK(  ) 

P 

logical  mask  for  nonzero  row  and  column  coordinates  in 
WORK(  ) 

LQWK(  ) 

P 

LIFO  queue  of  nonzero  coordinates  in  WORK!  ) 

PO(  ) 

P 

next  basic  factored  row  in  pre-order 

P(  ) 

P 

off-diagonal  row  with  nonzero  factored  coefficient 

D(  ) 

P 

depth,  remaining  back-substitution  path  length  in  factored 
component 

GN 

VGN(  ) 

P 

generalized  network  cycle  factors 

JMUL(  ) 

n 

ratio  of  generalized  network  coefficients  in  each  column 

Generalized  upper  bound  (GUB).pure  network  (PN),  and  generalized  network  (GN)  factorizations  respectively 
require  more  data  structures  to  support  kernel  factorization.  For  each  factorization.  F22  is  maintained  in  some 
partial  ordering  of  rows,  and  of  columns  —  a  signed  identity  for  GUB.  upper-triangular  for  PN  (e.g,.  Bradley, 
Brown,  and  Graves.  1977),  and  nearly-upper-triangular  for  GN  (e.g.,  Brown  and  McBride.  1984)  Direct  solution 
of  factored  kernel  equations  F,:z2-h:  and  -h\  is  performed  alternately  using  the  data  structures  shown. 
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are  reported  where  diversity  of  size,  structure,  and  taxonomy  have  proven  interesting  to  the 
modelers.  For  those  models  that  employ  nonlinear,  mixed  integer,  or  decomposition  features, 
we  report  solution  statistics  for  the  initial  linear  program. 

•  GTE.  The  seven  Telephone  Operating  Companies  within  GTE  have  adopted  an  inte¬ 
grated  business  system  called  Capital  Program  Management  System  ( CPMS )  to  guide  their 
3  billion  dollar  per  year  capital  planning.  The  system  includes  a  large  scale  mixed  integer 
programming  optimization  system  that  optimizes  the  critical  economic  tradeoffs  between 
maximizing  the  long-term  budget  value  of  the  firm’s  equity  and  satisfying  shorter-term 
financial  constraints,  resource  limitations  and  service  objectives.  Investment  opportunities 
for  the  next  5  years  are  modeled  as  0-1  variables  with  alternative  implementations  for  each. 
The  objective  is  to  maximize  the  net  present  value  of  the  capital  portfolio.  There  are  financial 
constraints  on  capital,  internally  generated  funds,  net  income  to  common,  and  limits  on 
resources  such  as  labor  hours,  lines  installed,  etc.  There  are  also  constraints  that  enforce 
logical  relationships  among  opportunities  (such  as,  if  choose  A  then  must  choose  B).  See 
Bradley  ( 1986). 

•  INVEST.  Capital  allocation  and  project  selection  for  Mobil  Oil  Corporation  are  modeled 
as  a  two-stage  multi-year  nonlinear  capital  budgeting  problem  with  over  40  000  integer 
variables.  A  master  problem  allocates  capital  among  markets  over  a  multi-year  horizon 
considering  the  estimated  nonlinear  effects  on  sales  of  concentrated  marketing  investments. 
The  instance  reported  here  is  a  mixed-integer  linear  program  subproblem  of  the  two-stage 
model  which,  given  these  annual  capital  expenditure  limits  for  a  market,  selects  particular 
alternate  investments.  Such  subproblems  are  easy  to  solve,  and  optimality  is  achieved  with 
a  single  iteration  of  the  nonlinear  master  problem.  See  Harrison,  Bradley  and  Brown  ( 1 989 ) . 

•  TANKER.  A  crude  oil  tanker  scheduling  problem  faced  by  a  major  oil  company  is 
solved  using  an  elastic  set  partitioning  model.  The  model  takes  into  account  all  fleet  cost 
components,  including  the  opportunity  cost  of  ship  time,  port  and  canal  charges,  demurrage 
and  bunker  fuel.  The  model  determines  optimal  speeds  for  the  ships  and  the  best  routing  of 
ballast  (empty)  legs,  as  well  as  which  cargoes  to  load  on  controlled  ships  and  which  to  spot 
charter.  All  feasible  schedules  are  generated,  the  cost  of  each  is  accurately  determined  and 
the  best  set  of  schedules  is  selected.  See  Brown,  Graves  and  Ronen  ( 1987). 

•  HFDF.  A  large-scale  elastic  set  partitioning  model  used  to  assign  frequencies  for  a 
network  of  high  frequency  direction  finding  receivers.  See  Brown,  Drake,  Marsh,  and 
Washburn  ( 1990). 

•  GAS.  A  multi-time  period  strategic  model  for  use  by  natural  gas  utilities  for  determining 
optimal  contract  levels  for  gas  purchase,  storage  and  transmission.  An  underlying  general¬ 
ized  network  flow  model  represents  gas  being  bought,  stored,  shipped  and  consumed  over 
a  multi-year  horizon,  typically  at  a  monthly  level  of  detail.  Constraints  and  variables  are 
added  to  handle  variable  maximum  and  minimum  purchase  levels,  variable  leased  or  con¬ 
structed  storage  and  variable  transmission  capacities.  An  integrated  parallel  model  incor¬ 
porates  the  peak  requirements  necessary  on  some  days  during  cold  winter  months.  This 
model  has  been  used  by  a  number  of  utilities  including  Southwest  Gas  Corporation  and 
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Questar  Pipeline  Corporation  to  plan  operations  and  to  justify  such  plans  to  regulatory 
agencies.  See  Avery,  Brown,  Rosenkranz  and  Wood  ( 1992). 

•  KELLOGG.  A  multi-time  period,  multi-plant  production/inventory/  transshipment  lin¬ 
ear  program  for  Kellogg  cereals.  The  model  guides  weekly  processing,  packaging  and 
shipping  decisions.  Production  consists  of  two  stages:  processing  lines  produce  basic  prod¬ 
ucts  which  are  then  packaged  on  packaging  lines  into  different-sized  containers  to  yield 
finished  products.  Processing  lines  produce  a  subset  of  the  basic  products  and  have  limited 
capacity  with  overtime  charges  for  weekend  shifts.  Packaging  lines  are  analogous.  In-house 
inventory  capacity  is  limited  although  outside  storage  is  available  to  additional  cost.  Inter¬ 
plant  shipments  of  finished  products  are  made  by  rail  or  truck.  See  Wood  ( 1989). 

•  ODS.  A  commonly  occurring  problem  in  distribution  system  design  is  the  optimal 
location  of  intermediate  distribution  facilities  between  plants  and  customers.  A  multicom¬ 
modity  capacitated  single-period  version  of  this  problem  is  formulated  as  a  mixed  integer 
linear  program.  A  solution  technique  based  on  Benders  Decomposition  is  developed — An 
essentially  optimal  solution  is  found  and  proven  with  a  surprisingly  small  number  of  Benders 
cuts.  See  Geoffrion  and  Graves  (1974).  The  instances  reported  here  are  decomposition 
master  problems. 

•  TAM .  The  annual  decision  on  how  much  the  Air  Force  should  spend  on  aircraft  and  on 
munitions  is  of  great  interest  to  many  people.  How  the  Air  Force  staff  develops  information 
to  support  the  decision  has  changed  over  the  years.  Currently,  a  linear  program  is  being 
used  by  the  Air  Force  Center  for  Studies  and  Analysis  and  is  being  tested  by  the  Munitions 
Division  of  the  Plans  and  Operations  Directorate  (AF/XOXFM)  for  munitions  tradeoff 
analysis.  The  LP  uses  existing  data  and  estimates  of  ( 1 )  aircraft  and  munition  effectiveness, 
(2)  target  value,  (3)  attrition,  (4)  aircraft  and  munition  costs,  and  (5)  existing  inventories 
of  aircraft  and  munitions.  Other  factors  considered  are  weather  and  length  of  the  conflict. 
See  Might  ( 1987)  and  Jackson  ( 1989). 

•  PHOENIX.  A  planning  model  for  the  multi-year,  multi-billion  dollar  modernization  of 
the  U.S.  Army’s  aging  helicopter  fleet.  The  mixed  integer  linear  program  employs  a  multi¬ 
product  production/inventory  formulation  with  aged  inventory.  Goal  constraints  attempt  to 
enforce  fleet  size,  maximum  age,  and  technology  goals  for  each  year  and  each  of  four  aircraft 
missions,  while  also  keeping  expenditures  within  upper  and  lower  limits.  Additionally, 
combinatorial  constraints  and  variables  handle  production  line  startup  and  shutdown  costs, 
minimum  and  maximum  production  levels  and  requirements  linking  certain  production 
lines.  See  Clemence,  Teufert,  Brown  and  Wood  ( 1988)  and  Brown,  Clemence,  Teufert  and 
Wood  (1990). 

•  EA6B.  Configures  jammers  of  hostile  radar  on  an  EA-6B  “Prowler”  Naval  electronic 
warfare  aircraft.  See  Sterling  ( 1990). 

•  DEC.  Digital  Equipment  Corporation  uses  this  model  to  determine  worldwide  manu¬ 
facturing  and  distribution  strategy  for  new  products.  This  mixed  integer,  linear  program 
suggests  a  production,  distribution,  and  vendor  network  which  minimizes  cost  and/or 
cumulative  cycle  times  subject  to  constraints  on  estimated  demand,  local  content,  and  joint 
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capacity,  over  multiple  products,  echelons,  and  time  periods.  Cost  factors  include  fixed  and 
variable  production  charges,  distribution  via  multiple  modes,  taxes,  duties  and  duty  draw¬ 
back,  and  inventory  charges.  See  Harrison,  Amtzen  and  Brown  ( 1992). 

•  AMMO  4H.  A  four-commodity  transshipment  model  for  delivery  over  time  of  military 
products  from  production  and  storage  locations  to  overseas  locations  to  support  theater 
operations  is  developed.  The  model  covers  five  physical  echelons,  including  production 
plants,  storage  depots,  ports  of  embarkation,  ports  of  debarkation  and  geographic  field 
locations.  Road,  rail,  sea  and  air  transportation  are  modeled,  and  product  demands  are  time- 
phased.  Capacitation  occurs  primarily  on  sea  and  air  links,  and  as  throughput  capacities  on 
transfer  points,  requiring  replication  of  some  echelons.  The  objective  of  the  model  is  to 
minimize  deviation  from  on-time  deliveries.  See  Staniec  ( 1984) . 

•  BUSCH.  A  model  of  brewery-to-wholesaler  movements  of  beer  for  Anheuser-Busch. 
The  model  also  includes  some  packaging  decisions  and  is  essentially  a  multicommodity 
flow  model  with  joint  capacity  constraints  arising  from  loading  dock  and  inventory  capac¬ 
ities  as  well  as  some  managerial  requirements.  See  Brown,  Mamer,  McBride  and  Wood 
( 1992).  The  instance  reported  here  is  a  small  pilot  model  for  the  full-scale  system  with 
millions  of  variables  which  is  solved  directly,  or  by  decomposition. 

•  BAR.  A  linear,  mixed-integer  multi-period  production-inventory  master  planning 
model.  See  Harrison  ( 1992). 

Four  implementations  are  compared:  “XS”  is  an  unadorned  version  of  the  X-System, 
an  implementation  of  the  Graves  mutual  primal-dual  method  with  its  GUB  factorization 
disabled,  while  “XS(GUB)”,  “XS(PN)”,  and  “XS(GN)”  each  employ  the  respective 
factorizations  discussed  here.  To  establish  a  frame  of  reference,  performance  of  these 
implementations  is  compared  with  two  well-known  commercial  solvers:  IBM’s  Optimiza¬ 
tion  Subroutine  Library  “OSL”  (Release  2,  1991 ),  and  two  versions  of  “CPLEX”  (Ver¬ 
sion  1.2,  1990,  and  Version  2.0,  1992). 

Ideally,  one  would  develop  four  equivalent  formulations  of  each  model,  each  customized 
for  its  particular  solver  with  the  goal  of  inducing  a  large  factored  row  set  of  the  appropriate 
type.  This  approach  is  a  consistent  theme  in  the  literature  dealing  with  specialized  algorithms 
and  one  that  we  strongly  endorse.  Alternate  formulations  of  a  model  are  often  available, 
and  it  seems  sensible  to  choose  one  that  exploits  as  much  as  possible  the  strengths  of  the 
solver. 

However,  all  of  the  models  used  here  are  “off-the-shelf”  in  the  sense  that  they  were 
developed  at  various  times  by  various  modelers,  and  alternate  formulations  are  impracti¬ 
cable.  Thus,  the  approach  is  to  preserve  a  single,  unfactored  representation  of  each  model, 
and  attempt  to  identify  favorable  row  structures  through  the  use  of  heuristics.  The  procedure 
is  based  on  the  work  of  Brown  and  Thomen  ( 1 980) ,  Brown  and  Wright  ( 1983 )  and  Brown, 
McBride  and  Wood  ( 1985).  The  heuristics  are  greedy  and  myopic  in  the  sense  that  they 
initially  consider  the  entire  row  set  of  the  problem,  and  discard  one  row  at  a  time  without 
backtracking  until  the  remaining  set  satisfies  the  desired  row  factorization.  This  can  be 
expected  to  confound  or  destroy  structure  introduced  by  the  modeler.  Although  the  auto- 
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matic  factorization  implementation  has  options  to  accept  modeler  guidance,  the  methods 
are  compared  here  without  this  subjective  complication.  While  these  model-naive  experi¬ 
ments  yield  interesting  and  useful  observations  about  the  implementations,  they  suffer  for 
lack  of  guidance  by  a  skilled  modeler. 

Table  9.1  shows  the  important  structural  information  concerning  the  model  instances  to 
be  solved. 

Table  9.2  displays  solution  times  for  the  sample  problems.  These  CPU  times  exclude 
initial  problem  input,  factored  row  identification,  and  final  output  —  on  average  about  0.2 
second  per  problem. 

The  original  formulations  of  most  of  the  test  problems  are  strongly  influenced  by  the 
modeler’s  solution  strategy.  For  instance,  TANKER  is  endowed  with  a  GUB  structure 
which  places  every  binary  variable  in  an  associated  set  from  which  only  one  member  can 
be  chosen;  this  maximal  GUB  set  is  also  sought  for  its  tendency  to  yield  nearly-integer 


Table  9.1 

Problem  dimensions 


71 

m 

Pgub^ 

Pgn^ 

NZEL 

GTE 

6  624 

960 

909/95 

909/95 

922/96 

58 

INVEST 

II  989 

1  338 

941/70 

1  101/82 

1  168/87 

33 

TANKER 

7  598 

83 

32/39 

32/39 

66/80 

31 

HFDF 

10  548 

61 

31/51 

31/51 

32/52 

189 

GAS  PN  A 

27  884 

6  848 

4  345/63 

5  934/87 

5  976/87 

37 

GAS  PN  C 

15  362 

3  794 

2  658/70 

3  418/90 

3  420/90 

20 

GAS  PN  E 

5  102 

1  184 

434/37 

877/74 

883/75 

7 

GAS  GN  A 

27  884 

6  848 

4  484/65 

5  142/75 

5  976/87 

37 

GAS GN  C 

15  362 

3  794 

2  664/70 

3  084/81 

3  420/90 

20 

KELLOGG  2 

17  841 

3  818 

1  265/33 

2  578/68 

2  596/68 

35 

KELLOGG  3 

27  490 

5  727 

2  295/40 

3  867/68 

3  892/68 

54 

KELLOGG  4 

37  139 

7  636 

2  428/32 

5  156/68 

5  188/68 

74 

KELLOGG  5 

46  788 

9  545 

3  388/36 

6  445/68 

6  484/68 

93 

ODS  1 

1 1  568 

3  023 

528/17 

540/18 

558/18 

21 

ODS  3 

23  993 

594 

490/82 

490/82 

490/82 

68 

TAM  5 

10  531 

438 

102/23 

132/30 

162/37 

94 

TAM  8 

6  104 

420 

118/28 

154/37 

196/47 

49 

TAM  12 

17  793 

629 

177/28 

231/37 

294/47 

165 

PHOENIX  10 

6  884 

1  618 

206/13 

220/14 

1  153/71 

14 

PHOENIX  30 

17212 

4  305 

293/07 

303/07 

3  604/84 

48 

EA6B 

12  247 

2  978 

1  610/54 

2  921/98 

2  921/98 

16 

DEC 

14518 

2  171 

677/31 

677/31 

1  088/50 

24 

AMMO  4H 

83  497 

13  963 

6  874/49 

12  892/92 

12  892/92 

129 

BUSCH  4 

7  997 

1  248 

649/52 

1  140/91 

1  148/92 

15 

BAR 

49  032 

7  446 

2  712/36 

4  575/61 

5  134/69 

102 

For  each  problem,  the  total  number  of  structural  variables  is  n.  structural  constraints  m.  GUB  rows  found  by  the 
identihcation  heuristic  pure  network  rows  pPS.  generalized  network  rows  pGS,  and  thousands  of  nonzero 
technological  coefficients  NZEL.  For  example.  GTE  has  58  thousand  nonzero  technological  coefficients  for  6  624 
variables:  GTE  can  be  viewed  as  having  960  explicit  and  no  factored  rows,  or  as  909/95%  GUB-factored  and 
960  —  909  =  5 1  explicit  rows,  as  909  PN-factored  rows,  or  as  922  GN-factored  and  38  explicit  rows. 
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Table  9.2 
Solution  seconds 


AMDAHL  5995-700  486/33  MHz  PC 


X-System 

IBM 

XS 

CPLEX 

none  GUB 

PN 

GN 

OSL 

GN 

1.2 

2.0 

GTE 

9 

8 

8 

8 

43 

41 

65 

58 

INVEST 

4 

5 

4 

4 

23 

24 

52 

49 

TANKER 

8 

10 

.10 

8 

6 

48 

8 

9 

HFDF 

100 

101 

101 

111 

40 

462 

581 

542 

GAS  PN  A 

2  376 

778 

50 

57 

291 

321 

1  714 

1  221 

GAS  PN  C 

4 

4 

3 

3 

79 

26 

185 

245 

GAS  PN  E 

72 

33 

4 

6 

13 

33 

165 

50 

GAS  GN  A 

1  115 

542 

294 

72 

312 

385 

3  532 

2  262 

GAS  GN  C 

4 

4 

3 

3 

71 

26 

186 

236 

KELLOGG  2 

3 

2 

2 

2 

69 

5 

219 

194 

KELLOGG  3 

5 

5 

5 

5 

144 

9 

513 

440 

KELLOGG  4 

48 

36 

26 

24 

500 

115 

1  274 

1  ill 

KELLOGG  5 

I  122 

1459 

320 

248 

1  210 

926 

2  615 

2  243 

ODS  1 

6 

3 

2 

5 

125 

22 

1  042 

414 

ODS  3 

6 

6 

6 

6 

23 

38 

58 

56 

TAM  5 

55 

60 

45 

40 

44 

231 

151 

86 

TAM  8 

24 

14 

14 

17 

21 

69 

77 

71 

TAM  12 

108 

112 

88 

106 

101 

312 

416 

378 

PHOENIX  10 

4 

2 

2 

2 

20 

10 

84 

73 

PHOENIX  30 

41 

25 

26 

9 

335 

69 

1  171 

1  239 

EA6B 

75 

33 

9 

9 

45 

44 

177 

244 

DEC 

189 

32 

42 

80 

71 

93 

317 

293 

AMMO  4H 

'  42 

41 

35 

46 

1  000 

277 

1  762 

1  597 

BUSCH  4 

5 

5 

4 

4 

26 

34 

55 

46 

BAR 

290 

212 

106 

114 

382 

480 

1  366 

1222 

An  AMDAHL  5995-700  running  under  IBM  VM/CMS/XA  with  IBM  VS  FORTRAN  2.3.0  is  used  to  render 
performance  in  CPU-seconds  accurate  to  the  precision  shown  for  the  basic  “XS”-system  using  no  factorization, 
compared  with  dynamic  factorizations  of  “XS(GUB)”,  pure  network  "XSf PN)”,  and  generalized  network 
“XS(  GN)”  rows.  “IBM-OSL”  shows  primal  simplex  performance  on  the  same  computer  of  the  IBM  Optimi¬ 
zation  Subroutine  Library,  Release  2  ( 1991 ).  “486/33”  shows  the  clock-time  performance  of  “XS(GN)”  on  a 
microcomputer  ( 33  MHz  Intel  486  with  32  MB  RAM)  followed  by  that  of  “CPLEX”  (Version  1.2)  (1990)  and 
of  “CPLEX”  (Version  2.0)  ( 1992)  on  the  same  microcomputer. 


solutions  to  the  linear  program.  AMMO  4H  is  a  multicommodity  capacitated  transshipment 
problem  and  thus  is  best  suited  to  a  pure  network  factorization;  it  was  originally  solved  by 
dual  decomposition  rendering  pure  network  sub-problems.  PHOENIX  is  a  multicommodity 
equipment  replacement  model  closely  following  the  generalized  network  factorization  par¬ 
adigm. 

The  row  structures  in  Table  9.2  have  been  found  without  modeler  help  by  our  automatic 
identification  heuristic,  yet  they  corroborate  the  modelers’  intentions,  TANKER  reveals  the 
same  pure  network  rows  as  the  GUB  rows,  or  a  generalized  network  constructed  by  iden- 
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tifying  one  additional  row  to  be  paired  with  each  GUB  row.  PHOENIX  exhibits  a  dominant 
structure  that  is  clearly  a  generalized  network. 

One  would  anticipate  the  factorization  exploiting  the  dominant  row  structure  to  win 
computation  tests.  This  is  wrong  more  often  than  right.  Table  9.2  shows  that  the  more 
general  factorization  dominates  the  less  general,  with  few  exceptions:  GTE’ s  relatively  large 
GUB  set,  and  the  large  pure  networks  in  GAS  PN  A,  GAS  PN  E,  and  AMMO  4H  seem  to 
satisfy  our  prior  bias  toward  model-dictated  factorizations. 

Table  9.2  also  suggests  that  our  myopic  use  of  heuristics  to  automatically  identify  factored 
structure  has  its  pitfalls.  In  a  number  of  problem  instances,  we  identify  significantly  larger 
factored  sets  with  the  more  general  factorizations,  yet  we  enjoy  little  improvement  in 
computation  times  (e.g.,  INVEST,  ODS  3,  TAM  8).  This  suggests  that  the  “quality”  of  a 
row  factorization  is  not  completely  specified  by  its  size. 

We  have  pursued  this  notion  by  inviting  some  of  the  modelers  to  guide  our  identification 
heuristic  to  precisely  the  row  sets  they  intended.  Some  of  the  results  have  been  striking: 
Wood  ( 1989)  reports  significant  improvements  for  problems  in  the  GAS  system. 


Table  9.3 

Maximum  number  of  elements  in  explicit  transformation  kernel 


XS 

XS(GUB) 

XS(PN) 

XS(GN) 

GTE 

13 

0 

0 

0 

INVEST 

9 

1 

1 

1 

TANKER 

1 

1 

1 

0 

HFDF 

3 

1 

1 

1 

GAS  PN  A 

1  550 

891 

30 

54 

GAS  PN  C 

18 

5 

0 

0 

GAS  PN  E 

263 

110 

7 

5 

GAS  GN  A 

1  470 

758 

340 

59 

GAS  GN  C 

19 

7 

1 

0 

KELLOGG  2 

6 

2 

0 

0 

KELLOGG  3 

9 

4 

0 

0 

KELLOGG  4 

79 

41 

10 

11 

KELLOGG  5 

1  299 

1  015 

138 

128 

ODS  1 

9 

1 

1 

5 

ODS  3 

0 

0 

0 

0 

TAM  5 

28 

22 

15 

18 

TAM  8 

20 

15 

8 

10 

TAM  12 

38 

42 

16 

21 

PHOENIX  10 

32 

21 

21 

1 

PHOENIX  30 

169 

157 

185 

1 

EA6B 

1  155 

270 

0 

0 

DEC 

218 

39 

49 

46 

AMMO  4H 

235 

92 

0 

0 

BUSCH  4 

10 

5 

0 

0 

BAR 

290 

163 

70 

41 

The  number  of  nonzero  elements  ( in  nearest  thousands)  in  the  explicit  transformation  kernel  A  I  j ,  gives  some 
indication  of  how  much  information  is  not  captured  by  factorization,  as  well  as  an  idea  of  relative  storage 
requirements. 


G.G.  Brown,  M.P.  Olson  /  Mathematical  Programming  64  (1994)  17-51 


45 


Table  9.4 

Binding  explicit  constraints  at  optimality 


Binding/% 

GUB/% 

PN/% 

GN/% 

GTE 

554/58 

20/02 

20/02 

18/02 

INVEST 

763/57 

201/15 

195/15 

162/12 

TANKER 

51/61 

60/72 

39/47 

9/11 

HFDF 

58/95 

29/48 

29/48 

28/46 

GAS  PN  A 

3  068/45 

1  953/29 

299/04 

349/05 

GAS  PN  C 

2  324/61 

850/22 

68/02 

90/02 

GAS  PN  E 

901/76 

573/48 

90/08 

79/07 

GAS  GN  A 

3  064/45 

1  854/27 

1  155/17 

359/05 

GAS  GN  C 

2  323/61 

861/23 

402/11 

89/02 

KELLOGG  2 

1  950/51 

1  015/27 

92/02 

118/03 

KELLOGG  3 

2  942/51 

1  368/24 

148/03 

183/03 

KELLOGG  4 

4  136/54 

2  491/33 

391/05 

452/06 

KELLOGG  5 

5  270/55 

2  305/24 

592/06 

617/06 

ODS  1 

361/12 

83/03 

76/03 

117/04 

ODS  3 

410/69 

0/00 

0/00 

0/00 

TAM  5 

264/60 

227/52 

168/38 

186/42 

TAM  8 

279/66 

240/57 

142/34 

175/42 

TAM  12 

412/66 

352/53 

205/33 

251/40 

PHOENIX  10 

1  098/68 

1  096/68 

1  090/67 

80/05 

PHOENIX  30 

3  298/77 

3  236/75 

3  239/75 

110/03 

EA6B 

2  939/99 

1  334/45 

26/01 

27/01 

DEC 

I  328/61 

736/34 

726/33 

421/19 

AMM0  4H 

6  686/46 

3  153/22 

13/00 

14/00 

BUSCH  4 

840/67 

481/39 

5/00 

8/01 

BAR 

4  687/63 

3  250/44 

1  895/25 

1  450/19 

Not  all  constraints  are  binding  at  optimality.  The  first  column  lists  the  number  of  binding  explicit  constraints  and 
expresses  this  as  a  percentage  of  all  constraints;  the  following  columns  display  the  number  of  binding  explicit 
constraints  and  their  corresponding  percentages  under  the  alternate  factorizations. 


A  few  of  our  corresponding  modelers  have  had  the  opportunity  to  build  models  from 
scratch  with  a  particular  factorization  in  mind.  This  admits  model  coercion  and  a  wide  range 
of  well-known  reformulation  methods  which  we  think  can  materially  change  both  the  size 
and  quality  of  the  result.  Their  early  reports  show  promise.  Among  the  models  discussed 
here,  the  GAS  and  KELLOGG  systems  have  been  subsequently  reformulated  to  enhance 
generalized  networks,  GTE  has  been  re-engineered  to  further  accentuate  its  dominant  GUB 
set,  and  TANKER-like  and  many  ODS  models  have  been  moved  to  a  micro-computer;  all 
these  models  are  now  larger,  but  much  easier  to  solve. 

It  is  surprising  and  encouraging  that  the  transition  to  more  general  factorizations  seldom 
degrades  performance  much,  even  when  few  additional  factored  rows  are  won  by  the 
increased  generality.  This  contradicts  popular  folklore  that  the  more  general  factorizations 
demand  substantial,  if  not  overwhelming  increases  in  the  resulting  sizes  of  the  factored 
structures.  In  fact,  computational  testing  reported  by  others  has  usually  been  limited  to 
models  in  which  the  number  of  explicit  rows  is  in  the  range  of  one  to  twenty  (e.g.,  Chen 
and  Saigal,  1977;  Glover.  Kamey,  Klingman  and  Russell,  1978;  Glover  and  Klingman, 
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1981).  Our  results  are  all  the  more  remarkable  given  the  lack  of  guidance  from  the  modeler 
for  the  “intended'’  row  factorization. 

Table  9.3  shows  the  maximum  size  of  the  A  f, 1  in  terms  of  its  nonzero  elements.  We  see 
that  the  maximum  size  of  the  explicit  transformation  kernel  tends  to  decrease  as  the  gen¬ 
erality  of  the  factorization  increases.  Recalling  the  definition  of  the  explicit  transformation 
kernel, 

A  1 1 1  =  ( E\ |  —  E\2F  22  ^21 )  '» 

this  trend  is  as  we  would  expect.  Each  potentially  binding  explicit  row  which  can  be 
converted  to  a  factored  row  reduces  the  likely  size  of  A  f, 1 ,  Also,  the  density  of  the  term 
-  El2  F  22 '  F2 1  generally  increases  with  the  density  of  F  22  •  For  F22  k-by-k,  the  number  of 
nonzeros  in  F  22  for  the  GUB  factorization  is  k,  for  pure  networks  perhaps  as  large  5  k2, 
and  for  generalized  networks  as  large  as  k2.  There  are  some  exceptions  to  this  trend  in  Table 
9.3,  especially  in  model  instances  in  which  the  size  of  the  (GN)  factored  row  set  is  not 
significantly  larger  than  that  of  (PN).  This  is  because  for  a  given  factored  kernel  F22,  the 
exact  (PN)  representation  of  F 22  is  generally  more  sparse  than  that  of  the  less  exact 
floating-point  representation  of  ( GN ) . 

It  is  usually  the  case  that  many  constraints  are  not  binding  at  optimality,  as  can  be  seen 
in  Table  9.4.  A  distinguishing  feature  of  dynamic  factorization  is  the  ability  to  limit  attention 
to  binding  constraints,  handling  binding  factored  constraints  with  great  efficiency,  and 
working  with  a  relatively  small  number  of  binding  explicit  constraints.  Explicit  binding 
constraints  on  the  order  of  a  few  thousand,  or  less,  are  quite  manageable.  This  is  well  beyond 
the  size  of  previously  reported  implementations. 


10.  Conclusions 

Previous  research  by  others  generally  suggests  that  specialized  algorithms  such  as  those 
presented  here  are  useful  only  when  the  factored  structure  completely  dominates.  There  are 
even  reports  of  algorithms  for  solving  problems  having  a  single  unfactored  (explicit) 
constraint  (Hultz  and  Klingman,  1978;  Klingman  and  Russell,  1978).  When  implementa¬ 
tions  have  been  reported,  problem  suites  have  been  limited  to  instances  having  a  very  small 
number  of  explicit  constraints,  typically  in  the  range  from  one  to  twenty  (Chen  and  Saigal, 
1977;  Glover,  Kamey,  Klingman  and  Russell,  1978;  Glover  and  Klingman,  1981).  The 
consensus  seems  to  be  that  such  algorithms  are  quite  delicate,  and  deserve  to  be  viewed  as 
specialized  algorithms,  useful  only  for  solving  very  special  problem  instances. 

We  refute  this  view.  Dynamic  factorization  is  competitive  with  commercial-quality 
optimization  systems  on  every  model  instance  we  have  tested. 

The  development  here  stresses  the  similarities  among  the  algorithms  and  the  natural 
extensions  leading  from  one  to  the  next.  This  is  in  contrast  to  the  development  reported  for 
simitar,  non-dynamic  algorithms  (e.g.,  Dantzig  and  Van  Slyke,  1967;  Klingman  and  Rus- 
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sell,  1978;  and  Hultz  and  Klingman,  1978)  in  which  the  specifics  of  the  individual  algorithm 
obscure  the  generality  of  the  approach.  The  conceptual  difference  between  our  algorithms 
is  seen  to  be  largely  isolated  to  the  structure  of  a  single  algebraic  entity,  the  factored  kernel. 
By  abstracting  the  structure  of  the  factored  kernel  and  concentrating  on  the  general  algorithm 
design,  the  versatility  and  flexibility  of  this  approach  is  clarified. 

The  algorithmic  development  leads  directly  to  an  implementation.  The  resulting  software 
suite  exhibits  a  “single  system  image”.  The  modularity  of  the  algorithm  allows  the  defi¬ 
nition  of  an  “abstract  data  type”  (see,  e.g.,  Aho,  Hopcroft  and  Ullman,  1974)  which 
isolates  the  data  structures  and  update  procedures  for  the  factored  kernel  from  the  rest  of 
the  implementation.  Each  factorization  is  seamlessly  integrated  within  the  system  design. 

The  early  1 980s  produced  a  great  deal  of  research  on  automatic  identification  of  special 
structure  in  LP  models  ( see.  e.g.,  Gunawardane  and  Schrage,  1977;  Glover,  1980;  Schrage, 
1981;  Brown,  McBride  and  Wood,  1985;  and  Bixby  and  Fourer,  1986).  We  have  incor¬ 
porated  the  most  useful  of  these  ideas  into  our  implementation,  and  we  have  what  we  believe 
to  be  the  first  complete  implementation  which  supports  automatic  identification  of  factored 
row  sets.  This  capability  may  be  used  to  identify  new  factored  structure  or  to  validate  or 
augment  a  modeler-provided  recommendation.  When  faced  with  the  choice  of  either  solving 
an  unfactored  model  instance  or  automatically  identifying  a  factored  structure  and  then 
using  the  corresponding  solver,  our  results  show  that  the  latter  is  nearly  always  to  be 
preferred.  Modelers  have  conducted  extensive  additional  computational  experimentation 
with  the  X-System  not  reported  in  this  paper.  These  results  suggest  that  in  addition  to  the 
quantity  of  factored  rows,  the  quality  of  these  rows  influences  the  performance  of  factori¬ 
zation  algorithms.  While  not  well  understood,  it  is  clear  that  the  myopic  approach  of  our 
heuristics  is  no  substitute  for  the  modeler’s  guidance  in  identifying  factored  structure. 

Processing  networks  ( Koene,  1982)  are  network  models  which  allow  proportional  flow 
restrictions  on  the  arcs  entering  or  leaving  some  nodes.  One  formulation  of  such  a  model 
results  in  a  pure  or  generalized  network  structure  with  a  set  of  complicating  columns.  Chen 
and  Engquist  (1986)  propose  a  primal  partitioning  algorithm  for  solving  processing  network 
problems.  An  alternate  formulation  yields  a  pure  or  generalized  network  structure  with 
complicating  rows;  this  is  precisely  the  structure  dealt  with  here. 

The  multicommodity  capacitated  transshipment  problem  (MCTP)  has  been  the  subject 
of  much  research  over  the  years,  and  a  number  of  specialized  algorithms  have  been  proposed 
to  solve  it  (see,  e.g.,  Assad,  1978,  or  Kennington,  1978) .  The  usual  MCTP  formulation  is 
a  pure  network  which  each  commodity  uses  independently  in  its  own  flow  model,  but  with 
side  constraints  on  the  total  common  flow  of  all  commodities  over  some  of  the  network 
arcs.  The  side  constraints  form  a  GUB  row  set,  while  the  rest  of  MCTP  forms  a  pure 
network;  either  view  might  be  preferred  depending  upon  size  of  the  common  network,  the 
number  of  side  constraints,  and  the  number  of  commodities.  In  our  experience,  the  network 
factorization  usually  dominates  the  GUB  factorization,  and  the  pure  network  factorization 
presented  here  is  a  powerful  technique  for  solving  MCTPs.  As  an  experiment,  we  customized 
our  (PN)  implementation  for  MCTP  to  exploit  the  special  structure  of  the  explicit  side 
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constraints.  This  highly-specialized  implementation  performed  no  better  on  AMMO  4H, 
and  we  now  believe  that  this  would  be  true  for  most  MCTPs. 

There  are  problems  which  would  exhibit  a  large  factored  row  structure  if  not  for  a  set  of 
complicating  columns  (e.g..  see  Brown,  McBride  and  Wood,  1985).  One  would  expect  the 
structure  of  the  factored  kernel  to  be  dominated  by  that  of  the  predominant  row  structure, 
with  only  occasional  complications  due  to  the  exceptional  columns.  One  might  allow  for 
this  exceptional  structure  in  the  factored  kernel  by  identifying  it  “on-the-fly”  as  the  algo¬ 
rithm  progresses,  and  preserving  the  sanctity  of  the  core  factorization.  Though  conceptually 
simple,  some  iterations  of  this  algorithm  would  border  on  the  spectacular.  This  approach 
may  be  thought  of  as  a  hybrid  between  the  dynamic  factorization  developed  here  and 
dynamic  basis  triangulation  methods  (see,  e.g.,  Hellerman  and  Rarick,  1971,  1972;  Saun¬ 
ders,  1976,  and  McBride,  1980). 

Dynamic  extrinsic  factorization  is  subsumed  by  the  algorithms  presented  in  this  paper  if 
we  activate  functions  in  the  update  analogous  to  the  secondary  exchanges  now  employed. 
Essentially  all  that  has  to  be  done  is  ensure  that  successive  factored  components  retain  their 
stipulated  special  structure.  We  speculate  that  this  will  work  best  in  cases  where  model 
structure  is  amenable,  and  quite  likely  will  require  some  model-specific  customization  to 
perform  well  on  difficult  models.  We  have  limited  our  experimentation  to  those  static 
extrinsic  cases  which  we  believe  to  be  most  generally  useful. 
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